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O ABSTRACT 

Context. The close-in planet orbiting GJ436 presents a puzzling orbital eccentricity (e = 0.14) considering its very short orbital 
period. Given the age of the system, this planet should have been tidally circularized a long time ago. Many attempts to explain this 
were proposed in recent years, either involving abnormally weak tides, or the perturbing action of a distant companion. 
Aims. In this paper, we address the latter issue based on Kozai migration. We propose that GJ 436b was formerly located further away 
from the star and that it underwent a migration induced by a massive, inclined perturber via Kozai mechanism. In this context, the 
perturbations by the companion trigger high amplitude variations to GJ 436b that cause tides to act at periastron. Then the orbit tidally 
shrinks to reach its present day location. 

Methods. We numerically integrate the 3-body system including tides and General Relativity correction. We use a modified symplec- 
tic integrator as well a fully averaged integrator. The former is slower but accurate to any order in semi-major axis ratio, while the 
latter is first truncated to some order (4 th ) in semi-major axis ratio before averaging. 

Results. We first show that starting from the present-day location of GJ 436b inevitably leads to damping the Kozai oscillations and 
to rapidly circularizing the planet. Conversely, starting from 5-10 times further away allows the onset of Kozai cycles. The tides act 
in peak eccentricity phases and reduce the semi-major axis of the planet. The net result is a two fold evolution, characterized by two 
phases: a first one with Kozai cycles and a slowly shrinking semi-major axis, and a second one once the planet gets out of the Kozai 
resonance characterized by a more rapid decrease. The timescale of this process appears in most cases much longer than the standard 
circularization time of the planet by a factor larger than 50. 

Conclusions. This model can provide a solution to the eccentricity paradox of GJ 436b. Depending on the various orbital configura- 
tions (mass and location o the perturber, mutual inclination. . . ), it can take several Gyrs to GJ436b to achieve a full orbital decrease 
and circularization. According to this scenario, we could be witnessing today the second phase of the scenario where the semi-major 
axis is already reduced while the eccentricity is still significant. We then explore the parameter space and derive in which conditions 

■ this model can be realistic given the age of the system. This yields constraints on the characteristics of the putative companion. 

Key words. Planetary systems - Methods: numerical - Celestial mechanics - Stars: Gliese 436 - Planets and satellites: dynamical 
C^l 1 evolution and stability - Planet-star interactions 
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1. Introduction balance leads for to an equilibrium temperature for the planet 

T e q — 642 K, assuming T e $ = 3350 K for G J 436 and zero albedo 

■ The M-dwarf GJ 436 has been the subject of growing interest foi ? ^ lanet Accordi to lDeming et al] rf2007h . the tempera- 
in recent years This st ar is known to ho st a close-in Neptune- ^ difference indicates tidal headi of GJ 436b but this could 
mass planet (Gl 436b Butler et alj | 2004) that was furthermore alternativel be due t0 the 8 ling the atmosphere in a 
observed to undergo transit d Gillon et alJ ] 2007|) . hot band if ^ lanet atmospher e does not radiate like 



X 



The monitoring of the transits of GJ 436b helped constrain- black bodv 
. ing its orbital solution. Noticeably, it appears to have a sig - 



nificant non-zero eccentricity: e = 0.14 ± 0.1 (Maness 2007 
furthermore refined to e — 0.14 ± 0.01 (iDemorv et alJ 12007.. . 
This eccentricity is abnormall y high for a small period planet Man y theones were proposed in recent years to account for 
(P - 2.64 days Ballarf!3 £010). With such a small or- the resldual eccentricity of GJ436b. The most straightforward 
bital period, tidal forces are expected to circularize the orbit one 18 * at the tldes m sufficiently weak and/or the age of the 
within much less than the present age of the system (1-10 Gyr s y stem ls sma11 enou 8 h t0 allow a re 8 ular tldal circularization 
Torres 2007). Tidal forces seem indeed to be at work in this n °U° be achieved yet. Thisjr^ea was_proposed by | Mardling| 
system. The detection of the secondary transit of the planet <™> for GJ436b, after |Trjlling| <|2000|) had suggested that 
enabled Demorv et al. (2007) to derive a brightness tempera- thls accounts for any close-in exoplanet observed with signa- 
ture of T = 712 ± 36 K, with flux reradiated across the day- lcant eccentricity. It is in fact related to our poor knowledge 
side only. Conversely, a stellar irradiation / thermal re-radiation of the P lanet ' s tldal dissipation Q„. Q p is a dimensionless pa- 

rameter related to the rate of energy dissipated per orbital pe- 

Send offprint requests to: H. Beust riod by tidal forced oscillations (Bar ker & Qgilviell2009l) . The 

Correspondence to: Herve.Beust@obs.ujf-grenoble.fr smaller Q p , the more efficient the tidal dissipation is. To lowest 
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order, the circularization time-scale f C jrc of a n exoplanet reads 
dGoldreich & Soterill966Hjackson et alJl2008h 



,13/2 



tr.\rc. — 



63 ^GM, 



■ Q P m p R 5 p 



(1) 



where M» is the stellar mass, a is the planet's orbital semi-major 
axis, m p its mass and R p its radius. Assuming Q p = 10 5 , R = 
27,600km, a = 0.0287 AU, m = 23.2 M m and M = 0.452 M 
dMardlingj2 008). this formula gives f c ; rc = 4.7 x 10 7 yr, which is 
obviously less than the age of the system. But Q p is very badly 
constraine d. For Neptune, it is estimated between ~ 10 4 an d 
3.3 x 10 5 dBanfield & Murray! [1991 IZhang & Harrriltonl l2008). 
We can thus conside r 10 5 as a standard likely value for GJ 436b, 
but iMardlingl (2008) argues that Q p could be as high as a few 
10 6 . In this context, assuming the lower bound (1 Gyr) for the 
age of the system, the circularization could not be achieved yet. 
This depends h owever on the starting eccentricity at time zero. 
IWisdonl 12008) performed analytical calculations of energy dis- 
sipation rates for synchronized bodies at arbitrary eccentricities 
and obliquities. His concluded that for a large enough starting 
eccentricity 0.4, as will be the case below), the circular- 
ization time-scale should be significantly reduc ed with r espec t 
to Eq. ([TJ. In this context, the conclusions by Mardling (2008) 
may no longer hold. This shows also that conclusions drawn on 
asymptotic formulas like Eq. (Q]) must be treated with caution. 
Numerical work is required. 

Alternatively, many authors proposed that the eccentricity of 
GJ 436b is sustained by perturbations by an outer massive planet, 
despite a standard Q p value. Such a lo ng period companion was 
initially suspected as iManessI d2007l) 's radial velocity data re- 
vealed a linear drift. More data di d not confirm the trend, which 
is now believed to be spurious. Montagnie r et al.l d2012l) gives 
instead observational uppers limits for the putative companion 
that rule out any companion in the Jupiter mass range up to a 
few AUs. 

Perturbers can be either resonant or non-resonant. Looking 
for a for a planet locked in mean-motion resonance with GJ 436b 
appears a promising idea, as resonant configuration usually 
trigger larg e r ecce ntricity modulations. This was suggested by 
Ribas et al.l (2008), who fitted in the residuals of the radial ve- 
locity data a ~ 5M e planet locked in 2:1 mean motion reso- 
nance with GJ 436b. However, this additional planet was invali- 
dated furthermore. Dynamical calculations by Mardling (2008) 
showed that this planet cannot sustain the eccentric ity of GJ 436b 
unless Q p ^ 10 6 . Moreover. lAlonso et all (120081) showed that 
a ~ 5M e planet in 2:1 resonance with GJ436b should trig- 
ger transit time variations that should have been detected yet. 
Such variations are un seen today with a high level of confi- 
dence dPont et al. 2009). Actually the constraints deduced from 
radial velocity moni t oring, transit time and g eometry monitoring 
dBallard et al.l l2010: Mo ntagnier et al.ll2012h almost rule out an 
additional large planet locked in a low-order mean-motion res- 
onance with GJ436b (P < 8.5 days), and to have a Jovian like 
planet up to a f ew AUs. 

Conversely, Ton g & Zhoul (120091) investigated the effects of 
secular planetary perturbations in eccentricity pumping, consid- 
ering both non-resonant (secular) and resonant configurations. 
They found some perturbers configurations that could account 
for the present eccentricity of GJ436b while still being com- 
patible with the observational constraints, but with no tides. 
Incorporating the tides leads inevitably to damp all the ec- 
centricity modulatio n and to circularize the orbit. Meanwhile, 
Batv gin et al.l d2009l) reanalyzed carefully this issue. They found 



that as expected, in most cases the eccentricity of GJ436b 
quickly falls to zero thanks to tidal friction, but that this effect 
can be considerably delayed if the 2 planet system lies initially 
in a specific configuration where the eccentricities of the two 
planets are locked at stationary points in the secular dynamics 
diagram. In this case, they show that the circularization time can 
be as high as ~ 8 Gyr des pite a standard < 9 D valu e. 

To date the model by iBatvgin et al.l d2009l) is the only one 
able to explain the high present day eccentricity of GJ436b 
with a standard Q p value. This model could appear as not 
generic, as it requir e s a s pecific initial configuration, but 
Batvgin & Morbidelli (2011) showed that in the framework 
of Hamiltonian planetary dynamics with additional dissipative 
forces (which is the case here), these points tend to behave as 
attractors. Hence various initial configurations can lead to reach 
su ch a point with a fur ther dynamical evolution like described 
by [Batvgin et al. (2009). As is shown below, in the model pre- 
sented here, we exactly encounter such a configuration where 
many different routes can lead to such a stationary point (the 
transition between Phase 1 and Phase 2, see Sect. 4). 

The purpose of this paper is to present an alternate model 
based on Kozai mechanism, assuming again a distant per- 
turber. This model can be viewed as complementary to the 
Batv gin et al.l (120091) study. Kozai mechanism is a major dynam- 
ical effect in non-coplanar systems that can trigger eccentricity 
modulations up to very high values. We make a review on this 
effect, without and with tidal effects taken into account and de- 
scribe our approach in Sects. 2 and 3 respectively. In Sect. 4, we 
present an application to the case of GJ436. We first present cal- 
culations starting from the present day orbit of GJ 436b, with the 
result that even Kozai mechanism cannot overcome the damping 
effect of tidal friction. Then we present other simulations where 
GJ 436b is initially put much further away from the star than to- 
day. We show that Kozai mechanism pumps GJ 436b's eccentric 
regularly up to high values where tides become active at perias- 
tron. The result is a decay of the orbit that drives it to its present 
day location but on a time scale that can be considerably longer 
than the standard circularization time. Hence even after several 
Gyrs, GJ436b's eccentricity can still be significant. In Sect. 5, 
we describe a parametric study of this scenario and derive in 
which conditions it is compatible with the present day situation 
of GJ 436b. Our conclusions are presented in Sect. 6. 

2. Kozai mechanism 

2.1. General features 



Kozai mechanism was first described bv lKozail (1 1962b . initially 
for comets on inclined orbit (S 50°) with respect to the eclip- 
tic. Under the effect of secular planetary perturbations (mainly 
arising from Jupiter), their orbit is subject to a periodic evolu- 
tion that drives it to lower inclination but very high eccentricity, 
while the semi-maj or axis remains constant. As pointed out by 
Bail ev etaU dl992). this mechanism is responsible for the origin 
of most sun-grazing comets in our Solar System. In the circular 
restricted 3-body problem, the Kozai Hamiltonian that describes 
this secular evolution is obtained by a double averaging of the 
original Hamil tonian over the orbital m otions of the planet and 
of the particle dKinoshita & Nakaill 19991) . 

Assuming zero eccentricity for the primaries, the secular mo- 
tion of the particle is characterized, to quadrupolar approxima- 
tion (expansion of the secular Hamiltonian to 2 nd order in powers 
of the semi-major axis ratio) by the following constants of mo- 
tions dKinoshita & Nakail 1 1 9991: iKrvmolowski & Mazehl f 1 999: 
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Fig. 1. Contour plots of the Kozai Hamiltonian (|4j in a (a), e) 
space for H = 0.05 



iFord et alJ2000HBeust & Dutrevll2006l) 

F = (2 + 3e 2 )(3 cos 2 / - 1) + 15e 2 sin 2 /cos2w 
H = (1 -e 2 )cos 2 / 



(2) 
(3) 



where e stands for orbital the eccentricity of the particle, i for its 
inclination with respect to the orbital plane of the primaries and 
a) for the particle's argument of periastron relative to this plane. 
The first constant (F) is the reduced Hamiltonian itself and H is 
the reduced vertical component of the angular momentum which 
is a secular constant thanks to the axial symmetry of the averaged 
problem. Combining both we can eliminate the inclination 



F = (2 + 3e 2 ) 



3H 



-1+15* 



cos 2a> 



(4) 



Contour plots of this expression in (a>, e) space for a given H 
value are called Kozai diagrams. The Kozai mechanism regime 
is characterized by small values of H (high initial inclinations). 
Figure Q] shows such a diagram corresponding to H — 0.05. We 
note the two islands of libration around a> = 90° and u = 270° 
(or -90°) characteristic for Kozai resonance. In a pure 3-body 
Newtonian dynamics, the particle just secularly moves along one 
of these curves. It is then obvious to see that its evolution is char- 
acterized by strong periodic eccentricity increases. In the same 
time e increases, i decreases to keep H constant. If for a given 
particle a> librates around ±90°, it is said to be trapped in Kozai 
resonance. But even if at circulates, its evolution is character- 
ized by the same eccentricity / inclination modulation described 
above, the maximum e values being reached when co = ±90°. 

As quoted above, this regime is valid for small H values. If 
the p article has initial zero ecce ntricity, it can be shown analyti- 
cally (Ki noshita & Nakail 1999b that the minimum inclination re- 
quired for Kozai mechanism to start is arccos( y/3/5) - 39.23°. 
Alternatively, it can operate at smaller inclination if the initial 
eccentricity is larger. But a some point in the secular evolution, 
it will reach an inclination above this threshold. 

This dynamical m echanism i s also active in hierar- 
chical triple systems dHarringtonl 11968b ISoderhielml 19821: 
Eggl eton et all 119981: iKrvmolowsk i & Mazehl ll999i iFord et al l 
2000), the inner orbit being subject to a similar evolution un- 
der secular perturbations by the outer body. This is particularly 
relevant for triple systems because they are often non-coplanar. 



Similarly, it is expected to be active on planetary orbits in bi- 
nary star systems, as soon as the orbit of the planets are signifi- 
cantly non coplanar with that of the binary (i ^ 40°, see above). 
This can affect the stability of the planetary orbits, especially if 
sever al planets are present (llnnanen et al.lll997l: lMalmberg et al.l 
120071) . It was actually invoked to explain th e high eccentricity of 
some extrasolar planets in binary systems dHolman et alJI 1997c 
Libert & Tsiganisll2009l) . 



2.2. Kozai migration 

This mechanism could also apply to the GJ 436 case, as a way to 
sustain the eccentricity of GJ436b. But tides inevitably lead to 
a circularization within the standard timescale (see si mulations 
below and also application to the TyCra system in Beus t et al.l 
Il997h . 

To render Kozai mechanism active despite the presence of 
tides, one must lower the strength of the tides. This can be done 
assuming a higher Q p as suggested by Mardling (2008) (see 
Introduction), but this is even not enough. Even with an arbitrary 
large Q p , there is still forced apsidal precession due to the plan- 
etary and stellar quadrupole moments and to General Relativity 
(see Sect. 3). With the present day semi-major axis of GJ436b, 
this is enough to override the Kozai mechanism. 

Then only way to damp all these effects to allow the on- 
set of Kozai cycles is to assume a larger semi-major axis for 
the planet during a first evolutionary phase. The purpose of the 
present paper is to investigate this idea. Of course, the present 
orbital semi-major axis of GJ436b is well constrained from 
radial-velocity monitoring, but our idea is to suggest that it used 
to be larger in the past, typically 5-10 times larger than to- 
day. In this context, all perturbing effects listed are too weak 
to prevent the onset of Kozai cycles. Then when the eccen- 
tricity increases, the periastron drops, because as long as tides 
are not active, the semi-major axis remains a secular invari- 
ant. Hence tidal friction starts at periastron. This leads to a 
gradual decrease of the orbit (i.e., the semi-major axis) which 
causes it to shrink and to finally reach its present location. 
Th is process is called Kozai migration. It has been described 
bv lEggleton & Kiseleva-Eggletonl(l200ll):IWu & Murravl d2003l) 
and more recently bv lFabrvckv & Trem aine (120071) . It is viewed 
today as one of the various migration mechanism that can act in 
planetary systems. Our aim is to investigate this scenario in the 
context of GJ 436b. We show that although it cannot prevent the 
final circularization of the planet, under suitable assumptions it 
can considerably slow down the circularization process, thus ex- 
plaining how the eccentricity of GJ 436b could still be significant 
a few Gyrs after the formation of the system in spite of strong 
tidal friction. 



3. Tidal forces and N-body integration 

3.1. Basic tidal forces 

Our goal is to test the scenario of Kozai migration to the GJ 436b 
case. This needs to be done via n umerical integration. T he vari- 
ous tidal forces are described in Mardling & Linl d2002l) . We re- 
call briefly this theory here. Let us consider i) a star with mass 
m s , radius R s , an apsidal constant k s and a dissipation parameter 
Q s , ii) a planet with mass m p , radius R p , an apsidal constant k p 
and a dissipation parameter Q p . The tidal forces acting on the 
planet are 
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1. the acceleration due to the quadrupole moment of the planet 
due to the distorsion produced by the presence of the star: 



7QD,p 



R 5 p (l + m s /m p )kp 



r 4 



5 (Sl p ■ uf 



\2Gm s 



-i(a p -u)si p \ , (5) 



where £l p is the spin vector of the planet, r is the distance 
between the planet and the star, and u — r/r, a unit vector 
parallel to the radius vector of the orbit; 
2. the acceleration produced by the tidal damping of the planet 



7tf, p 



6nk p m s I R„\ 5 /a\ 8 r n 



where a is the semi-major axis and n the mean motion of the 
orbit. 

Of course we have symmetric accelerations 7qd, s and 7tf,s 
for the star. We must add to these forces the relativistic post- 
newtonian acceleration 



7iel 



G (m p + m s ) 



(1 + 3rf)r ■ r - 2(2 + rf) 



G (m p + m s ) 



2 W 



u - 2(2 - rf)rr\ 



where rj = m p m s /(m p + m s ) 2 and c is the speed of light. Finally, 
the equation of motion describing the evolution of the orbit of 
the planet is 

G (m s + mp) 

r = ; r + y rel + y QD , p + 7 TF;P + y QDiS + 7 TFjS (8) 

r 

We must also add an equation go verning the evolution o f the spin 
vectors. For the planet we have dMardling & Linl l2002) 



m s m p i \ 

I p Sl p = r x 7qd,p + 7tf,p 

m s +m p 



(9) 



and a symmetric equation for Sl s . Note that we typically have 
Tqd.s « TQD.p and 7 TF , S « 7tf, p by a factor m p /m s , showing 
that the star itself is very little affected. To lowest order for a 
synchronized planet, we have 7qd,p/7tf,p — H2p/6. With Q p ^ 
10 5 . We see that 7qd, p dominates the orbital evolution. However, 
7tf,p (the damping term) is responsible for the circularization 
and must be retained. Also, we have to lowest order 



7rel 



r ■ r r m D 



7qd,p H k p c 2 Rl m s 



(10) 



With the numerica l values for the parameters of GJ 436b quoted 
in Mardling (2008), this ratio is 1.6. This shows that relativis- 
tic effects are expected to be comparabl e to tides and must not b e 
neglected. However, as pointed out by Mardlin g & Linl (|2002), 
the only secular effect of the relativistic potential is to produce 
an apsidal advance. This also holds for 7qd, p as both effects are 
conservative. They do not affect the circularization process due 
to 7tf, p - 



3.2. Inclusion into N-body integration and averaging 

As we want to investigate the combined effect of tides and plan- 
etary perturbations, these forces must be included into a N- 
body integration scheme. Our basic scenario is an inner planet 
(GJ 436b) subject to tidal interaction from the star, and perturbed 
by an outer, inclined and more massive planet (with mass m! and 
semi-major axis a') that is itself not subject to tides. 

A major difficulty is that we have to deal with phenomena 
with very different timescales. In the context of GJ436, the pe- 
riod of Kozai cycles is typically ~ 10 5 yr, while the tidal circu- 
larization acts on timescales ^ 10 8 yr. Hence an efficient A^-body 
integrator is required. 

A standard way to overcome this difficulty is to average all 
interaction terms over the orbital motions of both bodies. This 
can concern the interaction 3-body Hamiltonian between the 
planet and the perturber, as well as the tidal terms. This way 
only the secular part of the perturbations is retained. 

The averaging of the tides and of the relativistic terms is 
done only over the orbital period of the inner planet. To do 
this, it is convenient to derive the effect of the tides on the 
specific angul ar momentum vector h = r X r and the Runge- 
Lenz vecto r (Mardling & Lin 2002 cjFabrvckv & Tremainel 2007: 
Eggleton & Kiseleva-Eggleton 2001) 



rxh 



G (m s + nip) 



(ID 



(n\ Each additional acceleration y will induce a rate of change of e 
and h given by 



dh de 2 (7 • r)r - (r ■ r)y - (y ■ r)r 
— -rxy ; — = (12) 



dr 



dt 



G (m s + m p ) 



These equations are then averaged over the orbital period of 
the inner planet. Using Hansen coefficients, th is can be done 
in clo sed form. Explicit formulas are given by Mard ling & Linl 
(I2002h . 

The averaging of the interaction Hamiltonian with the outer 
perturber over the orbital period of both bodies is more complex, 
as this cannot be done in closed form. One has to expand the 
corresponding terms in ascending powers of the semi-major axis 
ratio a I a' between the planet and the perturber (using Legendre 
polynomials). Assuming a/a' <K 1, it is tr uncated to som e given 
ord er and then averaged oyer both orbits. IWu & Murravl (l2003h 
and iFabrvckv & Tremaind d2007l) truncate the Ha miltonian to 
second order (quadrup olar), while Ford et ail ( 20001) recommend 
third order (octupolar). Mardling &Linl(l2002h~ give explicit for- 
mulas up to third order (for coplanar orbits though). As ex- 
plained below, the level of accuracy required by our study forced 
us retain terms up to 4 th order in a /a', i.e. one level beyond oc- 
tupolar. We used Maple to derive analytic formulas up to this 
level. 

Another way to handle the 3-body interaction is to use 
symplectic integration. The main difficulty is to combine it 
with tidal forces. The theory of symplectic integration is based 
on the Hamilton equations of motion that apply to any N- 
body problem. I t s back grou nd is for i nstanc e described in 
ISaha & Tremaind dl992h and IChambersI dl999h . The key idea 
is to split the Hamiltonian H into pieces, say H — Ha + Hb, 
each of them being exactly integrable, and then to integrate each 
part s eparately alternatively. Most common sy mplectic integra- 
tors (fChambers Sl999l : lLevison & Duncan! 19941) are based on the 
following second order scheme: At each time-step r, we must 



4 



H. Beust et al.: Dynamics of the Gliese 436 system 



1. advance Hb by a half step t/2; 

2. advance Ha for the full time step r; 

3. Re-advance Hb by t/2. 

In planetary and hierarchical A^-body problems, Ha corresponds 
typically to a collection of independent Keplerian Hamiltonians 
and Hb to the disturbing function. We will name this scheme 
B(t/2)A(t)B(t/2). Our goal here is to include tides to it, i.e. 
dissipative forces. Of course with dissipation the integration, as 
well as the real problem, is no longer symplectic. This does not 
prevent the integration method from being efficient, at least for 
the A^-body forces. Typically, with the second order scheme de- 
scribed above, a time-step of ~ 1 /20 of the smallest orbital pe- 
riod is enough in standard problems t o ensure a c onservation 
of the energy to ~ 1CT 6 relative level dBeustll2003h . Including 
dissipative forces leads to a non preservation of the energy, but 
the same time-step can be kept. There have been actually many 
successful attempts in recent years to include dissipative forces 
into symplectic integrators (Qvlalhotra 19941 iKehoe et al.| [2003; 
ICuk & Burnl2004HMarzari & Weidenschilling||2006l) . 

The most straightforward way to proceed is to include the 
integration of the various forces given above into the B(t/2) sub- 
steps (the disturbing part of the integration). However this turns 
out not to be efficient in our case. The reason is that the tidal 
forces are highly dependant on the distance r between the planet 
and the star (oc r 4 and even steeper). Therefore, when consider- 
ing significantly eccentric orbits, keeping the standard symplec- 
tic time-step leads to an insufficient mapping of the tidal forces 
close to periastron. 

To overcome this difficulty, we used averaged formu- 
las only for the tidal and relativistic terms. Then we add 
to our symplectic scheme two sub-steps T(t/2) of integra- 
tion of these averaged equations for e and h, in the order 
B{t/2)T(t/2)A(t)T{t/2)B(t/2): 

1. advance Hb by a half step t/2; 

2. integrate the averaged tidal equations by t/2; 

3. advance Ha for the full time step r; 

4. re-integrate the averaged tidal equations by t/2; 

5. re-advance Hb by t/2. 

At each T(j/2) sub-step, we integrate simultaneously the evolu- 
tion of the spin vectors of the bodies via Eq. ((9]). The advantage 
of this approach is that it avoids the use of truncated expansions, 
while still being efficient. We modified this way o ur symplecti c 
integrator HJS dedicated to hierarchical systems dBeustl 120031) . 
All computations presented in Sect.|4]have been done using this 
integrator. Its use is nevertheless more computing time consum- 
ing than a fully averaged integrator. This is why we also built a 
fully averaged integrator for our parametric study (Sect. 5). 

To check the validity of our integration method, we also 
made some comparisons with integrations performed with a con- 
ventional Bulirsh & Stoer integrator. Of course the latter integra- 
tions were much more time consuming so that they were done 
over a much less extended time-span. In all cases the agreement 
appeared excellent. An example is shown on Fig.[3](eccentricity 
plot) where we have overplotted in red the result of the conven- 
tional integration. Both curves are almost identical, so that the 
red curve barely appears under the black one. 

4. Application to Gliese 436 

4.1. Integration starting from the present state 

The first element we want to test is how planetary perturbations 
may affect the present orbit of GJ436b. As mentioned above, 



given the present semi-major axis a = 0.0287 AU, tides are ex- 
pected to dominate the dynamical evolution. We therefore per- 
formed four runs starting from the present-day orbital configu- 
ration with different assumptions: 

- run (a): with a distant perturber, but not taking into account 
the tides nor the general relativistic post-newtonian correc- 
tion (hereafter GR), i.e. using the standard HJS code; 

- run (b): taking into account tides and GR, but with no per- 
turbing planet (run (b) in Fig. [2j; 

- run (c): with the same perturbing planet as in run (a), but 
taking now into account the relativistic post-newtonian cor- 
rection, but no tides; 

- run (d) : with the same planet taking into account both GR 
and tides. 

In all cases, the perturbing planet was chosen with mass m' - 
0.1 Mj up and semi-major axis a' = 3.583 AU (i.e., with orbital 
period T' = 10 yr) as a typical convenient set of parameters (see 
Sect. 5), compatible with the constraints of iMontagnie r et al.l 
(2012) based on radial velocity residuals and photometry. In ad- 
dition we fix the initial tilt angle ;'o between the two orbits to 75° 
to be able to generate a strong Kozai resonanc e. We assume the 
stellar, planetary and tidal parameters given by Mardling (2008) 
and listed in the Introduction. 

In Fig. [2] we show the secular evolution of the eccentricity 
of GJ 436b over 2 X 10 7 yr in the four cases. In case (a), we see 
as expected a large amplitude modulation characteristic for the 
Kozai resonance created by the perturbing planet. In case (b), 
we note a gradual circularization with a time-scale ^ 2 x 10 8 yr. 
In case (c) (same as (a) but with GR taken into account), there 
is no Kozai resonance anymore. The relativistic precession ac- 
tually smoothes the Kozai mechanism. The eccentricity is ba- 
sically constant with high frequency, small amplitude oscilla- 
tions. In run (d) (with tides and relativity), the high frequency 
oscillation is superimposed to a tidal decrease at the same rate 
as in case (b). The main conclusion we can derive from these 
runs is that the perturbing planet could in principle induce a 
strong Kozai resonance, but that tides and GR override this ef- 
fect. Even GR alone (run (c)) is enough to cancel out the Kozai 
resonance. The reason is that the secular precession induced by 
G R is larger than that indu ced by the Kozai resonance. As shown 
by [Fabrycky & Tremaind (120071) . GR adds an extra term to the 
secular precession of the argument of periastron co, and subse- 
quently, the maximum eccentricity of the Kozai cycles is sig- 
nificantly reduced. This can be understood considering Eq. (0). 
With an extra precession for a>, the cos 2a> term in the averaged 
Hamiltonian F can be considered as rapidly varying, so that it 
secularly averages to zero. Only the first term can be kept in F. 
This makes the Hamiltonian independent of all angles, making 
all related actions constant, and hence e and i are constant. This 
is what we get here. 

If we consider the dissipative case with tides, the net result in 
our case is that the tidal circularization process remains almost 
unchanged (run (d)) with respect to the case with no perturbing 
planet (run (b)). 

The Kozai resonance appears thus unable to sustain the high 
present eccentricity of GJ 436b over the age of the system. The 
only ways to solve this paradigm are either to increase the 
strength of the Kozai resonance either to lower the tides and the 
relativistic effects. The former way would imply to increase the 
mass of the perturbing planet by at least two orders of magni- 
tude; this would render it incompatible with the constraints of 
Montagni er et al.l d2012l) . The latter way can be achieved assum- 
ing a larger semi-major axis for GJ 436b. Of course the present 
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Fig. 2. Secular evolution of the eccentricity of GJ 436b starting from the present-day orbit and assuming a standard perturber (m' = 
0. lMj U p, T = 10 yr) in different conditions (see text): (a) perturbing planet, but no tides nor relativity; (b) tides and relativity but no 
perturbing planet; (c) perturbing planet, relativity but no tides, (d) perturbing planet, relativity and tides. 



semi-major axis is well constrained, but we may assume that it 
used to be larger in the past. This is why we come to the idea of 
Kozai migration. 

4.2. The Kozai migration scenario 
4.2.1. General description 

We describe now into details a typical run where we assume an 
initial semi-major axis for GJ436b significantly larger than the 
present value, and where both GR and tides are taken into ac- 
count, with unchanged tidal parameters with respect to Sect. 14. II 
We assume for the perturber the same parameters as in Sect. I4.ll 
(a' = 3.583 AU, m! = 0.1 M Jup , i = 75°), but we take as ini- 
tial semi-major axis for GJ 436b ao = 0.287 AU (i.e, 10 times its 
present-day value). 

Fig. [3] shows the secular evolution of GJ436b in this con- 
text over 2 x 10 7 yr. We see that the eccentricity is subject to a 
large amplitude modulation characteristic for Kozai resonance. 
With flo = 0.287 AU, apsidal precession arising from GR and 
tides is no longer able to override the Kozai resonance. But at 
peak eccentricity (e 0.8) in the Kozai cycles, the periastron 
drops to 0.05-0.06 AU, i.e., only twice the present-day semi- 
major axis. As a result, tides are at work in peak eccentricity 
phases. This can be seen in Fig. [3] We note a rapid reduction 
of the obliquity of GJ 436b relative to its orbital plane (initially 



fixed to 30°), so that after a few 10 7 yr, the spin axis of GJ 436b 
remains perpendicular to its orbital plane. We also note a spin- 
ning up of the planet's rotation. On Fig.[3](bottom-left), we plot 
as a function of time the ratio of the rotation velocity GJ 436b 
to its orbital mean-motion (Q/n), so that 1 means synchronism. 
This ratio was initially fixed to 2, and we note a gradual in- 
crease up to 10 at t = 2 x 10 7 yr. This reveals a spinning up of 
the planet, as the orbital mean-motion (i.e. the semi-major axis) 
is not subject to so drastic changes in the same time. Finally, 
we note a very small decrease of the semi-major axis over the 
time-span considered. These effects are clearly due to the tides 
at work in high eccentricity phases. This shows up obviously 
in the staircase-like shape of the temporal evolution of the pa- 
rameters displayed. Actually they are subject to some evolution 
only in the peak eccentricity phases, and remain unchanged in- 
between. The spinning up of the planet is also a consequence 
of the tides in high eccentricity phases. The rotation of GJ436b 
tries indeed to tidally synchronize with the orbital angular veloc- 
ity at periastron, which is very different from the mean-motion 
due to the high eccentricity. With a peak eccentricity of ~ 0.8, 
synchronism with angular velocity at periastron means actually 
Q/ra =! 16. Finally, the semi-major axis decrease reveals a loss 
of orbital energy due to tidal friction. 

We come now to describing the same simulation, but over a 
much longer time span, i.e. 20Gyrs (Fig. 0)). This timescale is 
intentionally taken larger than the age of the system and even of 
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Fig.3. Secular evolution over 2 X 10 7 yr of GJ436b assuming an initial value a = 0.287 AU and a perturber with m! = 0.1Mj up 
and a' = 3.583AJ7 (T' = 10 yr): Eccentricity (top left), semi-major axis (top right), rotation speed with respect to the orbital 
mean-motion (£2/n, bottom left), and obliquity angle (bottom right). On the eccentricity plot we overplot in red the result of the 
same integration, but using a conventional Bulirsh & Stoer integrator 



the Universe to show the possible future evolution. We display as 
in Fig.[3]the eccentricity, the semi-major axis, the Q/« ratio, and 
the mutual inclination between the two orbits. On the eccentric- 
ity plot, the color curves are calculated using various averaged 
integrators (see below). The one corresponding to our symplec- 
tic integration in the black one. We clearly see a two-fold evolu- 
tion, first before lOGyrs and then after. In the first phase, Kozai 
resonance is still active as can be seen from the high amplitude 
eccentricity oscillations. In the same time, the semi-major axis is 
subject to a slight decrease. The eccentricity oscillations do not 
remain identical though. Actually the peak eccentricity slightly 
increases but a more striking fact is that the bottom eccentric- 
ity of the Kozai cycles gradually increases to finally reach the 
peak eccentricity after ~ lOGyr. At this point (beginning of the 
second phase) the eccentricity oscillations stop, which means 
that the GJ436b gets out of the Kozai resonance. The semi- 
major axis drops then much more quickly and the eccentricity 
decreases to zero within a few extra Gyrs. During almost all the 
evolution (apart from the beginning), the rotation of GJ 436b re- 
mains synchronized with the angular velocity at periastron in 
high eccentricity phases. The increase of Q/n during the first 
phase is related to the gradual increase of the peak eccentricity. 
In the second phase O/n gradually decreases towards 1. At the 
end of the simulation, GJ 436b is tidally synchronized and cir- 



cularized, but the final semi-major axis is several times smaller 
than the initial one. The Kozai migration is ended. 

This two-f old behaviour in Koza i mi gration process was al- 
ready noted by Wu & M urravl (120031) and Fabrvc kv & Trem aine 
d2007l) . Actually ou r results o f Fig.|4]are v ery sim ilar to the cor- 
responding ones of iFabrvckv & Tremaind d2007l) . 

4.2.2. The first phase 

During the first phase, the planet still undergoes Kozai oscil- 
lations, but the bottom eccentricity of the cycles gradually in- 
creases to finally reach the peak one. This can be explained by a 
shift in Kozai diagram. 

In a pure 3-body Newtonian dynamics, the planet just sec- 
ularly moves along one of the level curves in a Kozai diagram 
similar to that of Fig.Q] In the w-circulating case, the curves are 
explored from left to right with dw/df > 0, and in the w-librating 
case (the genuine "resonant" orbits), the curves are explored 
clockwise around the stable point. Now, bo th GR and tides add 
an ex tra, positive component to dw/df dFabrvckv & Tremaind 
2007) at peak eccentricity in the cycles. This results in a slight 
shift towards right in the Kozai diagram each time the planet 
passes through a peak eccentricity phase. The planet jumps to 
another nearby curve located to the right. It is then obvious 
to see from Fig. Q] that this change induces a decrease of the 
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Fig. 4. Same evolution as in Fig. [3] but over 2 x 10 10 yr: Eccentricity (top left), semi-major axis (top right), rotation speed with 
respect to the orbital mean-motion (Q/n, bottom left), mutual inclination between the two orbits (bottom right). On all plots, the 
black curve corresponds to our symplectic integration. The color curves on the eccentricity plot are calculated with a fully averaged 
integrator, with interaction Hamiltonian truncated up to various orders: 2 nd order (blue), 3 rd order (green), and 4 th order (red). 



bottom eccentricity in the w-librating case and an increase in 
the w-circulating case. Now it is also obvious that thanks to 
these shifts in diagram at each eccentricity peak, a planet ini- 
tially in the w-librating case will be pushed to move to ca- 
circulating after a while. Consequently the bottom eccentric- 
ity of the planet is expected to eventually first decrease during 
the w-librating phase and then to further increase in the subse- 
quent ti»-circulating re g ime. T his is exactly what is reported by 
Fabry cky & Tremaine (2007) in their Fig. 1. After a short de- 
crease, the bottom eccentricity of the cycles increases gradually 
to 1. In our Fig. [4] we see that the bottom eccentricity (black 
curve) of the cycles does not immediatel y start to increase, but 
we do not see an initial decrease like in Fabrvckv & Tremaine 
J2007h . In fact, the more realistic integration is actuall y ours. The 
reason is that the integration of Fabr vckv & Tremaind d2007l) 
was made on the basis of a quadrupolar expansion of the 
Hamiltonian, as well as the Kozai diagram in Fig. Q] It is thus 
not surprising to see them in perfect agreement. But our integra- 
tion, made using a symplecti c integrator , is not truncated and is 
accurate to any order in a /q'. lFordet al. (2000) showed that re- 
taining higher order terms, in particular up to the octupolar level 
(3 rd order in expansion in a /a') considerably improves the accu- 
racy of the integration, especially in the vicinity of the separatrix 
between the two regimes, which is the case discussed here. 



In any case, when the planet enters the w-circulating regime, 
the bottom eccentricity of the cycles is forced to increase. 



4.2.3. The se cond phase and the link with the lBatvain et al.l 
12009) model 

When the bottom eccentricity reaches the peak one, then the 
planet moves along a more or less straight line in the Kozai di- 
agram. Tides and GR are now active permanently as the eccen- 
tricity is permanently high. This causes the precession of ai to 
increase. This forces the planet to start a more sharp tidal de- 
crease that makes it get out of the Kozai resonance. In this sec- 
ond phase, the planet gets gradually tidally circularized. 

An important outcome of this scenario is that it can provide 
a solution to the eccentricity paradox of GJ436b. We note in 
Fig. |4] that the first phase of the secular evolution lasts several 
Gyrs. Hence GJ436b may have first spent a long time in the 
first phase before starting a sharp tidal decrease. Moreover, we 
note that in the second phase, even if the semi-major axis shrinks 
quickly, it still takes a few more Gyrs to the eccentricity to de- 
crease downto zero. This could appear surprising, as we are now 
away from Kozai resonance. In Fig. |4] the situation of the planet 
at ~ 13 Gyrs is similar to the one we studied in Fig. [2] In the 
latter case, the eccentricity of GJ 436b decreases downto zero in 
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~ 2 x 10 8 yr, but in Fig.|4]it still takes a much longer time from 
13 Gyrs to reach zero. How could we explain this discrepancy ? 

The difference comes from the dynamical configuration 
of the 3-body system. It is in fact closely related to th e 
quasi-stationary sol ution described by Batvgi net alJ (|2009). 
Baty gin et al.l ([2009) studied the coplanar 3-body system. They 
showed that there exist two apsidal fixed points in this system, 
characterized by de/df = 0, and tu-vj' = orw-m' = n, where 
m and m' are the longitude of periastron of GJ436b and the 
perturber respectively. There is always one stable point among 
these two configurations around which librations can occur. In 
that case, m - m' just librates around its equilibrium value, and 
de/df = periodically vanis hes. As a result , the d ecrease of the 
eccentricity is much slower. Bat ygin et al.l (l2009h showed that 
adding tides to this peculiar system considerably increases the 
tidal circularization time, because the angular momentum trans- 
fer associated with tidal dissipation is now shared among the 
two interacting orbits. This is of course a peculiar situation, but 
as mentioned above, these points tend to behave as attractors. 

The situation we describe here differs from that of 
Baty gin et al.l (2009) because we do not have a coplanar config- 
uration. But the starting point of the second phase is also char- 
acterized by a stationary point with de/df = 0. At the end of the 
first phase (Kozai cycles), when the bottom eccentricity of the 
cycles equals the peak one, GJ 436b moves along a straight line 
in the Kozai diagram . Hence we have de/df = 0. Contrary to 
Baty gin et all <|2009), here m — vj' circulates during the second 
phase. The reason is that in the coplanar problem, there is no 
qu adruplar component to de/df. The expression of de/df given 
bv lBatvgin et al.l d2009b arises indeed from the octupolar expan- 
sion. But at high inclination, the quadrupolar component to de/df 
dominates, so that the constraints are different. But the common 
feature is that the starting point of both evolutions is character- 
ized by a stationary point in eccentricity in the pure secular 3- 
body evolution. The result is a huge increase of the tidal circular- 
iz ation time by a factor 50 at least, like in the evolution described 
by Batygin et al. (2009). Besides, the apsidal alignment criterion 
could be used as a diagnostic. If an apsidally aligned companion 
is discovered in the future in the radial velocity residuals of the 
star, it would be a strong in dication in favor of a coplanar con- 
figuration like described by Baty gin et alJ d2009l) : conversely if 
this is not the case, the companion could be considered as highly 
inclined with good probability. 



5. Exploration of the parameter space 

5. 1 . Description of the results 

We claim that the Kozai migration process we describe could 
fit the present day configuration of the GJ 436 system. GJ 436b 
could fairly well be today in the middle of the second phase, 
with an already reduced semi-major axis and a slowly decreasing 
eccentricity. In Fig. [4] the situation of the planet at ~ 13 Gyrs can 
be compared to the present day orbital configuration of GJ 436b. 

Of course 13 Gyrs is too long compared to the age of the 
GJ 436 system. But the timescale of the two-fold evolution can 
vary considerably depending on the initial parameters of the 
model, so that adjusting them to fit to any given age less than 
the age of the Universe is possible. 

We come now to exploring the parameter space of this sce- 
nario to try to derive in which case it can be compatible with the 
observation. We tested various orbital configurations, exploring 
the parameter space. The important parameters are 



- the mass m' and semi-major axis a' (or equivalently the or- 
bital period 7") of the perturber; 

- the initial semi-major axis ao of GJ 436b; 

- the initial mutual inclination j'o between the two orbits. 

Note that instead of z'o, the actual critical parameter is the re- 
duced vertical component of the specific angular momentum 
h = VT^ e 2 cos i, which is a secular constant of motion (See 
Sect. 2. 1). But in all our simulations we took an initial eccentric- 
ity eo = 0.01 for GJ 436b, so that h is unambiguously related to 

Other parameters such as the eccentricity of the perturber, 
the initial rotation speed and rotation axis orientation of GJ 436b 
turned out to have only a minor influence on the global behaviour 
of the system, so that exploring the space of the major parame- 
ters is enough for our purpose. In particular, in all the following, 
the initial eccentricity e' of the perturber was set to 0.03 as a 
typical standard value. 

The key outcome to monitor is the transition time between 
the two phases in the scenario (hereafter f tl ), which occurs close 
to 10 Gyrs in Fig. [4] This timescale appeared to be extremely 
sensitive to the parameters, ranging from less than 10 5 yr in some 
extreme cases to hundreds of Gyrs in other configurations. In all 
cases, most of the orbital decrease and circularization that occurs 
after f tr is achieved with ~ 1.5ftr- 

In our exploration, we consider as likely any t a value be- 
tween 1 Gyr and lOGyr. Shorter values are not compatible with 
a not yet circularized GJ 436b today, and with longer values, the 
star is not old enough to have entered the second phase yet. Our 
goal was to derive which sets of parameters yield a convenient t r . 
To do this, we first used our symplectic code in a few cases, and 
then our averaged integrator to save computing time. As men- 
tioned above, the averaging process of the Hamiltonian due to 
the perturber cannot be done in closed form. One has to first ex- 
pand it in ascending powers of a /a' , to tr uncate it to some given 
order and then average over both orbits. iFabrvckv & Trem aine 
(2007) truncate to second order (quadrupolar). while iFord et al. 
(2000) recommends third order (octupolar). We checked that be- 
ing in correct agreement with the symplectic integration (in par- 
ticular in the monitoring of f tl which is a very sensitive param- 
eter) requires to retain terms up to 4 th order in a/a', i.e. one 
level beyond octupolar. The use of a truncated integrator does 
not change the global two-fold behaviour reported above, but it 
affects its timescale. This is illustrated in Fig. [4^ (eccentricity 
plot) where we have superimposed to our symplectic integration 
(black curve) the same eccentricity evolution, but calculated with 
averaged integrators truncated at various orders: second (blue), 
third (green) and 4 th (red). The use of the quadrupolar trunca- 
tion generally leads to an error on f tl of ~ 50%. Most of the 
time it overestimates it. This can be understood easily: In a pure 
Newtonian dynamics (no tides, no relativity, such as in Fig.|2^), 
all Kozai cycles are strictly identical in the quadrupolar approx- 
imation. The motion is strictly periodic, as can be seen from 
Fig. Q] This is due to the absence of the perturber's eccentric- 
ity and argument of periastron in the interaction Hamiltonian 
(Eq. d2]>). This is no longer tr u e at higher order ; as shown by 
iKrvmolowski & Mazehl (119991) : IFord et al.l d2000l) . Due to sec- 
ular changes of the perturber's orbit, the Kozai cycles are no 
longer identical (see Fig. [3^ for instance). The total number of 
Kozai cycles achieved within a given (long) time-span can thus 
vary depending on the theory assumed. As the efficiency of the 
Kozai migration process depends on the integrated number of 
peak eccentricity phases reached, this explains the variability of 
f tl -. We see from Fig. |4^ that the use of a third order (octupolar) 
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theory (green curve) considerably enhances the agreement with 
the symplectic integration. A ~ 5% agreement is reached if we 
go to 4 th order (red curve). Or course better agreement could be 
reached if we used even higher order formulas. But in that case, 
the averaged formulas are so complex that the gain in computing 
time with respect to the symplectic integration is small. As a 5% 
accuracy on f tI is enough for our analysis, we decided to use 4 th 
order theory. 

Figure [5] shows a summary of the result of this exploration 
of the parameter space. Each plot corresponds to the use of a 
fixed set of values (ao, z'o) where we let the other parameters 
(T',m') vary. According to the simulation results, we divide the 
(T',m') plane into different regions depending on the resulting 
value of f tr . Regions marked "Fast Kozai" correspond to cases 
where t tt < 1 Gyr. In these regions, the Kozai migration pro- 
cess should have ended a long time ago, and we would expect 
GJ 436b to be synchronized and circularized today. This situa- 
tion logically applies to massive and/or close-in perturbers (top- 
left parts of the panels), when the disturbing action is stronger 
and the period of Kozai cycles smaller. Conversely, towards the 
bottom parts of the panels (small and/or distant perturbers), we 
have regions marked as "Slow Kozai", thus corresponding to 
f tr > 10 Gyr. In such cases, f tr is just too long to allow GJ436b 
to lie in the second phase of the Kozai migration process today. 
This cannot match the observation, as during the first phase of 
the evolution (Kozai cycles), the semi-major axis only undergoes 
a very small decrease an cannot reach the present day value. We 
also note in some cases regions marked "No Kozai" correspond- 
ing to very low mass and distant companions. In such cases, the 
perturbing action of the companion is so weak that GR (even at 
the considered distance) is strong enough to prevent the onset 
of Kozai cycles. Thus GJ 436b's eccentricity never grows and it 
can never reach its present-day location. To illustrate this fact, 
we add to all plots in Fig. [5] an oblique line which corresponds 
in each case to fcR = f Kozai, where ?gr and f Kozai are character- 
istic times for G R and for Kozai oscillat ions respectively. These 
timescales read (Matsumura et al. 2009) 
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Here we have computed these timescales assuming e — 0. Note 
that the real period of the Kozai cycles equals ZKozai within a fac- 
tor of order unity dFord et alJIiOOO: iBeust & Dutrevl l2006). It is 
straightforward to see that ?gr = ?Kozai implies m' oc T' 2 (with 
m' <s M), hence the oblique dotted lines in Fig. [5] No Kozai os- 
cillations are to be expected below this limit (?gr < z&>zai), which 
is confirmed by our numerical exploration. The "No Kozai" re- 
gions extend indeed somewhat above the dotted lines the our 
plots. This is due to the higher order effects that are not taken 
into account in the above simple formulas. 

In the middle of the fast and slow Kozai regions, a conve- 
nient strip appears where the value of ? tl allows GJ436b to lie 
presently in the middle of the second phase of the Kozai migra- 
tion process. Note that for a given set of parameters, the proba- 
bility of witnessing this phase today is not small, as the second 
phase lasts typically ~ 0.5?tr. 

To all plots in Fig. [5] we superimpose a s olid line that cor- 
respo nds to the observational upper limits of Montagnier et al. 
(2012) for any additional planet to GJ 436. The left oblique part 
corresponds to the radial velocity limit and the horizontal upper 
part to the imaging (adaptive optics) limit. Consequently, any set 



of (T', m') values located above this line must be rejected, and 
for an initial set of parameters (ao, z'o) to be possible, the conve- 
nient strip must cross the region below the observational limits. 
In all plots, the lower T value (in x-axis) was set according to 
the stability limit of the system. Systems with smaller T values 
are unstable. The stability was tested for these orbits with the 
unaveraged symplectic integrator. The upper limit in T was set 
at the point where the convenient strip definitely gets above the 
observational imaging limit, so that any convenient companion 
located further away would have been already detected. So, for 
any given set of initial conditions (ao, z'o), the convenient values 
for T and m! must appear on the plot. 

For all displayed cases, the convenient strip appears two- 
fold. Towards lower T values (close-in perturbers), the conve- 
nient m' values are a slowly decreasing function of T' , while for 
more distant perturbers, they increase rapidly. This can be under- 
stood with a careful analysis of the Kozai mechanism. The effi- 
ciency of the Kozai migration scenario depends actually on how 
close GJ 436b can get to the star during the first phase (hence the 
maximum eccentricity it can reach) and the time it spends at high 
eccentricity to allow tides to be at work. For close-in perturbers, 
the process is limited by the maximum eccentricity it can reach. 
In fact in Fig. [5] close-in perturbers in the convenient strip have 
masses comparable to that of GJ436b 0.074Mj up ), and have 
semi-major axes only a few times more than GJ 436b initially. 
Hence the initial angular momentum ratio I'J Iq between the per- 
turber and GJ436b is of order unity or a bit more. Now, due 
to total angular momentum conservation (we neglect exchanges 
with spins), we necessarily have 



l 2 Q + 1' + 2/q/q cos z'o = I 2 + I' 2 + 211' cos z 



(15) 



where I' and / stand for the angular momenta of both bodies at 
peak eccentricity, and i for the mutual inclination in the same 
phase. Assume now I' = l' Q (the perturber is not affected) and 

Z = l Vl -e 2 (the semi-major axis is only little affected in the 
first phase). This equation simplifies to 



2 j- (coszo • 



vr 



e- cos Z(> I + e 



) 







(16) 



Solving this equation for e shows that for moderate IUIq values, 
the eccentricity cannot go beyond a maximum value that is sig- 
nificantly below 1, while for larger value, any value as close pos- 
sible to 1 is allowed. In other words, the angular momentum of 
GJ 436b cannot reach zero because that of the perturber cannot 
grow to account for the whole angular momentum of the system. 
Let us assume for instance z'o = 75°, i = 30° (typical values) and 
Zq/Zo = 2. We derive e 0.74, which is high, but significantly be- 
low 1 to avoid GJ 436b to be efficiently tidally braked in peak ec- 
centricity phases. The efficiency of the Kozai migration process 
is therefore limited by the maximum eccentricity value than can 
be reached, hence by the value of Zq/Zo- Similar t a correspond to 
similar IUIq values. Assuming now l'„ oc m'T' l/3 and Zo oc mT 1 ? 3 
(low initial eccentricities), a constant f tl - leads m' oc T'^ 1 ^ 3 , hence 
the decrease of the convenient strip towards close-in perturbers. 

For more distant perturbers, the l' Q /lo ratio is large enough 
not to limit the peak eccentricity of the Kozai cycles. The ef- 
ficiency of the process is now limited by the amount of time 
spent in high eccentricity phases, as tides are only active there. 
We thus expect f tl - oc f KoZ ai, as the longer f tr , the less time is 
spent within a given time span in high eccentricity phases. As 
?Kozai K (Mjm!) x (T' 2 /T), a fixed value for f tl will correspond to 
a constant ratio T' 2 /m'. This explains the increase of the conve- 
nient strip towards distant perturbers. A more distant perturber 
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Fig. 5. Result of the exploration of the parameter space for four (flo,z'o) combinations, a: flo = 0.287 AU, z'o = 75°; b: ao = 
0.1435 AU, z = 75°; c: a = 0.287 AU, z = 85°; d: a = 0.574 AU, z'o = 85°. In each plot, the upper observational limits (radial 
velocities and imaging) are represented as a solid line; the oblique dotted line corresponds to the limit ?gr = fKozai as given by 
iMatsumura et al.l 12009). In the upper left plot, the dashed lined are m' oc T'-l/3 and to' oc T' 2 power laws (see text). 



just needs to be more massive to generate the same Kozai migra- 
tion as a closer one. 

To illustrate this simple analysis, we have drawn in Fig. [5^ 
two dashed lines corresponding to arbitrary m' oc j'- 1 ^ and 
m! oc T' 2 power laws. We note that these power laws described 
above as well followed in the simulations for both regimes, as the 
convenient strip quite faithfully remains parallel to the dashed 
lines sketched. To summarize, the convenient strip roughly fol- 
lows a dual power law regime, with m' oc j'- 1 ^ for small T 
and m' oc T' 2 for large T , Its exact location in (T',m') space 
depends on the initial set of parameters (ao, i'o), as can be seen 
from Fig. [5] 

We also note that in all cases, the convenient strip appears 
quite narrow, although this effect is enhanced on Fig. [5] by the 
logarithmic scale. In fact, for a given period T , Kozai migration 
is very sensitive to the mass to' of the perturber, so that to' must 
assume a well defined value within a factor ~ 1 .5 to allow f tI to 
fall in the right range. However, for a wide range of initial values 
sets (ao, z'o) (see discussion below), there is a convenient strip in 
(T', to') which is compatible with the observational constraints, 
so that the probability of the scenario is not small. 



5.2. Discussion 

We come now to compare the location of the convenient strip 
between various configurations. If we compare Figs.|5^ and [5J; 
(same ao but different z'o), we note that for z'o = 85° the conve- 
nient strip is located lower than for z'o = 75°. Here again this 
may be easily understood. With a larger z'o the Kozai mechanism 
is stronger, so that GJ 436b more easily gets a high peak eccen- 
tricity and the tides act more efficiently. Hence a less massive 
perturber is required to achieve the same result. 

If we compare configurations with same z'o but different ao 
(Figs. [5^ and b on the one hand, and c and d on the other hand) 
we can see how the convenient strip moves in the (T',m') dia- 
gram if we vary «o- We see that increasing ao brings the conve- 
nient strip up in the decreasing part of the diagram (m! oc r~ 1/3 ), 
and down in the increasing part of the diagram (to' oc T' 2 ). If for 
instance we superimpose Figs.|5^ and b, the convenient strip for 
ao = 0.1435 AU falls below that for a = 0.287 AU for small 7", 
then crosses it and get above for large T . This can be explained. 
In the decreasing part, the efficiency of the Kozai mechanism is 
controlled by the angular momentum ratio l' /k). for a larger ao, 
Iq is larger and the angular momentum ratio is smaller, so the 
Kozai mechanism is less efficient. We thus need a larger l' Q (and 
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hence a larger m') to get the same results. This is why the conve- 
nient strip is brought up for a larger ao- In the increasing part of 
the diagram, the Kozai mechanism is controlled by fKozai- As we 
know that fKozai K T 1 , a larger ao leads naturally to a smaller 
?Kozai and to a more efficient Kozai migration process. Hence a 
less massive companion (larger m') for the same T is required 
to achieve the same fKozai- This is why the convenient strip is 
brought down if we increase ao- 

This regime tends nevertheless not to be valid for all val- 
ues of ao- For larger initial orbits, the whole convenient strip 
is shifted upwards if we increase ao- For instance, for ao = 
0.574 AU and ;'o = 75° (configuration not displayed in Fig. [5]), 
the convenient strip is located much higher than with ao = 
0.287 AU, so that is it barely compatible with the observational 
constraints. Similarly, with io = 85°, the convenient strip is 
shifted upwards for ao ^ 0.6 AU. This behaviour is due to the 
fact that the efficiency of the tides are related to the minimum 
periastron reached rather than the maximum eccentricity. As in 
the first phase, the semi-major axis of GJ 436b is only little af- 
fected, for a larger ao a higher eccentricity is necessary to reach 
the required periastron range. 

Now, is it possible to set absolute limits to the unknown pa- 
rameters ? Setting a lower limit to ;'o is easy. As described above, 
if we decrease io the convenient strip gets higher in (T',m') 
space. At some point it crosses the observational limits so that 
the considered (ao, io) must be rejected. We tested several other 
runs with ;'o = 65° and io = 70°. In all cases the convenient 
strip (if any) appears above the observational limits. It is then 
possible to stress that io must be at least 70-75°. Note also that 
the Kozai mechanism is symmetrical with respect to io = 90°. 
We checked that runs assuming io = 95° and i'o = 105° behave 
exactly the same way as corresponding runs with io - 85° and 
io = 75° respectively. It is thus possible to put an upper limit 
to io to 105-110°. The initial configuration of the system must 
have been close to perpendicular within less than 20°. This is of 
course valid for small initial eccentricity eo- If <?o is high, then 
io may take different values. In the general case, a more relevant 
constraint must be derived for h = Vl — e 2 cos i. 70° ^ io 5 90° 
for small eo means \h\ S5 0.34, which is valid irrespective of eo. 

It is also possible to set a lower limit to ao. Looking at Fig. [5] 
we see that for lower ao the convenient strip gets shifted towards 
left, so that at some point it should cross the radial velocity ob- 
servational limit. This occurs for ao — 0.1 AU. Hence we can 
stress that the initial value ao must have been larger. 

Giving an upper limit to ao is less easy. We first can notice 
than choosing a larger ao implies a much more restricted compat- 
ible range for io. As noted above, for io = 75°, ao = 0.287 AU is 
possible, but ao = 0.574 AU is barely excluded (too high conve- 
nient strip). For io = 85° nevertheless, we tested that even fairly 
large ao values (up ~ 10 AU) we can find a convenient strip in 
(T',m') space that falls still below the observational limits. We 
may thus set an absolute upper limit to ao — 10 AU. However, we 
must point out that for such values, we are in an extreme Kozai 
regime with a very high peak eccentricity. Even if such a situa- 
tion cannot be rejected, whether the onset of such a regime with 
nearly orthogonal orbits is likely is questionable. Hence even if 
the absolute compatible range for ao is 0.1-10 AU, we may stress 
that the range 0.1-0.5 AU is more probable. 

As pointed out above, in all cases, the convenient strip in 
(T', m') space to allow f t , to fall in the suitable time range is quite 
narrow. This does not tell that the situation we are describing is 
improbable, as suitable configurations are possible for a wide 
range of initial configurations. This tells that for any given set of 



values (ao, io, T') the mass of the companion m' must fall around 
a well defined value within a factor ~ 1.5. We nevertheless note 
that taking larger ao values leads to generate narrower conve- 
nient strips. This shows up already in Fig.[5jl (ao = 0.574 AU), 
and it is even more striking for larger ao- This is another reason 
why we think that a configuration with a moderate ao is more 
likely. 

Our parametric study also shows that the evolution we de- 
scribe is not the only one possible. First the initial inclination 
io must be high to have a Kozai migration. Second, if the per- 
turber's mass m' is too small, there is no Kozai migration, and if 
it is too high, Kozai migration is extremely rapid. Kozai migra- 
tion was also suspected to account for the misaligned systems 
that have been measured from Rossiter-McLaughlin effect. But 
among all the systems for which an obliquity with respect to 
the stellar axis has been measured, only ~ 30 % show an obliq- 
uity larger than 45° (He brard et al.ll201 lb . This shows that Kozai 
migration is not a universal process. The limited size of the con- 
venient strip in our parametric study supports this idea. Systems 
like GJ 436 with eccentric close-in planets, although not unreal- 
istic, are not expected to be numerous. 

6. Conclusion 

The Kozai migration scenario is a generic process that causes 
inward migration of planets in non coplanar planetary systems. 
This evolution is characterized by two phases, a first one with 
perturbed Kozai cycles, and a second one with a more drastic 
decay of the semi-major axis. In any case, the characteristic time 
for this evolution can be considerably longer than the standard 
tidal circularization time, providing this way a solution to the is- 
sue of the high prese nt eccentricity of GJ 436 b. If we combine 
our study and that of Monta gnier et al.l d2012l) . various suitable 
initial configurations can be found. They all require a high initial 
mutual inclination (<: 75°), an initial orbit for GJ436b several 
times wider than today, perturber masses ranging between less 
than 0.1 Mj up and 50 Mj up , with orbital periods ranging between 
a few years up to several hundreds. This model is also generic 
enough to be invoked to explain some other puzzling cases of 
high eccentricity close-in exoplanets. But for a given set of ini- 
tial conditions, the orbital period and the mass of the perturber 
must assume a fairly well constrained relationship, otherwise the 
suspected process is either too rapid or nonexistent. This could 
then explain why system like GJ 436 are not numerous. 

Note that the constraint derived for io is valid for small ini- 
tial eccentricity eo. The real condition is \h\ $ 0.3, which can be 
achieved either with low eo and high ;'o, but also with low io and 
high eo. For small io, \h\ £ 0.3 implies eo ^ 0.95. Although this 
is actually just a matter of fixing the starting point on the Kozai 
cycles, this has a different meaning for what concerns the for- 
mation of such a planetary system. Planetary systems are usu- 
ally thought to form in circumstellar disk wi th initial coplanar 
configurations (see e.g. theoretical studies bv iPierens & Nelsonl 
2008: ICrida et alj 12009). The only way to imagine an initially 
inclined configuration is to assume that both planets formed in- 
dependently and were assembled further to build the GJ 436 sys- 
tem. This could happen for instance if GJ 436 was initially mem- 
ber of a multiple stellar system. Such systems are indeed rarely 
coplanar. Alternatively, the system could have formed copla- 
nar but the inner planet could have been driven to high eccen- 
tricity. There are many ways to achieve this. This can be done 
by pla net-planet scattering events or chaotic diffusion. Laskar 
(1 19961) indeed showed (with application to Mercury, see also 
Wu & Lifhwickl (1201 lh ) that in any planetary system, the inner 
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planet may be subject to drastic eccentricity increases due to sec- 
ular chaos. Similarly, recent work on the so-called Nice model 
for the formation of the Solar System dTsiganis et al.l 120051; 
iLevison et al.ll2.Ql T) have demonstrated the role of planet-planet 
scattering among the giant planets. Mean-motion resonances can 
also drive inner bodies to high eccentricities. This mechanism 
was invoked to explain the suspected Fallin g Evaporated Bodies 
(FEB) scenario in the /? Pictoris system (B eust & Morbidelhi 
ll996ll2Q00h . 

In conclusion, an initially nearly coplanar but eccentric con- 
figuration of the system appears more realistic than a highly in- 
clined one. Note that this does not change anything to the fur- 
ther evolution described in this paper, as Kozai cycles cause the 
system to rapidly oscillate between these two extreme configu- 
rations. Even if the system is initially nearly coplanar, it evolves 
to highly inclined within half a Kozai cycle. The relevant is- 
su e is actually the defin i tion o f the starting point. According 
to Batygin & Morbidelli (1201 ll) . in all systems potentially sub- 
ject to Kozai resonance, Kozai cycles are first frozen as long as 
self gravity and strong planet-planet interactions induce a large 
enough apsidal precession. But at some point when the system 
relaxes, Kozai cycles can start, which constitutes the "starting 
point" in our scenario. 

The mean-motion resonance model that applies to the FEB 
scenario can also provide a valuable starting point. The only re- 
quirement is that both planets are initially in a mean-motion res- 
onance like 4:1 ou 3:1, with a small initial mutual inclination 
(a few degrees) and a small eccentricity 0.05) for the outer, 
more massive planet. If this is fulfilled, the inner body is subject 
to an eccentricity incre ase virtually up to ~ 1 . It was shown by 
iBeust & Valironl (12007 ) that whenever it reaches the high eccen- 
tricity state (the FEB state in the case of /? Pictoris), it under- 
goes high amplitude inclination oscillations that can be viewed 
as Kozai cycles inside the mean-motion resonance. This is a kind 
of resonant version of our model that could constitute a full co- 
herent scenario driving the system from a coplanar and circular 
configuration to strong Kozai resonance. 

The reality of this scenario will nevertheless be difficult to 
test observationally. Even if the two planets are initially in mean- 
motion resonance, as soon the semi-major axis of GJ 436b starts 
to evolve thanks to tidal dissipation, they leave the resonance. 
The present day orbital configuration keeps no memory of the 
initial resonance. Assuming that we are today in the middle of 
the second phase with a still significant eccentricity and an al- 
ready reduced semi-major axis (by a typical factor 10), we are 
now unable to distinguish between initially resonant and non- 
resonant configurations. We also note from Fig.|4]that the present 
day mutual inclination between the two planets should be in any 
case ~ 20°, which makes i t difficult also to dis tinguish observa- 
tionally between coplanar dBatvgin et alJl2009h and Kozai mod- 
els (this work). However, as explained above, the exact char- 
acterization of the apsidal behavior of the planets can lead to 
constraints on their mutual inclination. Basically, if m - tzr' is 
close to or n, this would constitute a s trong clue in favor of the 
coplanar model of iBatygin et al.l {2009). Otherwize, this would 
indicate an inclined configuration. 

Before that, a crucial issue is to continue the radial velocity 
and photometric follow up of GJ 436. If we confirm the presence 
of a more distant companion, and if we are able to significantly 
constrain its orbit, then we will be able to test the proposed sce- 
narii into more details. This would not only resolve the eccentric- 
ity paradox of GJ 436b, but it would also provide clues towards 
the system's past evolutionary history. 



Acknowledgements. All computations presented in this paper were performed at 
the Service Commun de Calcul Intensif de FObservatoire de Grenoble (SCCI) 
on the super-computer funded by Agence Nationale pour la Recherche un- 
der contracts ANR-07-BLAN-0221, ANR-2010-JCJC-0504-01 and ANR-2010- 
JCJC-0501-01. We gratefully acknowledge financial support from the French 
Programme National de Planetologie (PNP) of CNRS/INSU. We thank our ref- 
eree, Konstantin Batygin, for valuable suggestions to improve the discussion of 
the paper. 



References 

Alonso R., Barbieri M., Rabus M., et al., 2008 

Bailey M.E., Chambers J.E., Hahn G., 1992, A&A 257, 315 

Ballard S., Christiansen J.L., Charbonneau D., et al., 2010, ApJ 716, 1047 

Banfield D., Murray N., 1992, Icarus 99, 390 

Barker A.J., Ogilvie G.I., 2009, MNRAS 395, 2268 

Batygin K., Laughlin G., Meschiari S., et al., 2009, ApJ 699, 23 

Batygin K., Morbidelli A., 2011, Celest. Mech. Ill, 219 

Beust H., Corporon P., Siess L., Forestini M., Lagrange A.-M., 1997, A&A 320, 
478 

Beust H., 2003, A&A 400, 1129 

Beust H., Dutrey A., 2006, A&A 446, 137 

Beust H., Morbidelli A., 1996, Icarus 120, 358 

Beust H., Morbidelli A., 2000, Icarus 143, 170 

Beust H., Valiron P., 2007, A&A 466, 201 

Butler R.P., Vogt S.S., Marcy GW, 2004, ApJ 617, 580 

Chambers J.E., 1999, MNRAS 304, 793 

Crida A., Baruteau C, Kley W., Masset E, 2009, A&A 502, 679 

Cuk M., Burns J.A., 2004, Icarus 167, 369 

Demory B.O., Gillon M., Barman T., et al., A&A 475, 1125 

Deming D., Harrington J., Laughlin G, et al., 2007, ApJ 667, L199 

Eggleton P.P., Kiseleva L.G., Hut P.W., 1998, ApJ 499, 853 

Eggleton P.P., Kiseleva-Eggleton L., 2001, ApJ 562, 1012 

Fabrycky D., Tremaine S., 2007, ApJ 669, 1298 

Ford E.B., Kozinsky B., Rasio F.A., 2000, ApJ 535, 385 

Gillon M., Pont F, Demory B.-O., et al., 2007, A&A 472, L13 

Goldreich P., Soter S., 1966, Icarus 5, 375 

Harrington R.S., 1968, AJ 73, 190 

Hebrard G, Ehrenreich D., Bouchy E, et al., 2011, A&A 527, Ll 1 

Holman M., Touma J., Tremaine S., 1997, Nature 386, 254 

Innanen K.A., Zheng J.Q., Mikkola S., Valtonen M.J., 1997, AJ 1 13, 1915 

Jackson B., Greenberg R., Barnes R„ 2008, ApJ 678, 1396 

Kehoe T.J.J., Murray CD., Porco C.C., 2003, AJ 126, 3108 

Kinoshita H., Nakai H., 1999, Celest. Mech. 75, 125 

Krymolowski Y., Mazeh T., 1999, MNRAS 304, 720 

Kozai Y, 1962, AJ 67, 591 

Laskar J., 1996, Celest. Mech. 64, 115 

Levison H.F., Duncan M.J., 1994, Icarus 108, 18 

Levison H.F., Morbidelli A., Tsiganis K., Nesvomy D., Gomes R., 201 1, AJ 142, 
152 

Libert A.-S., Tsiganis K., 2009, A&A 493, 677 
Malhotra R., 1994, Celest. Mech. 60, 373 

Malmberg D., Davies M.B., Chambers J.E., 2007, MNRAS 377, Ll 

Mardling R.A., Lin D.N.C., 2002, ApJ 573, 829 

Mardling R.A., 2008, arXiv0805.1928 

Marzari E, Weidenschilling S.J., 2006, BAAS 38, 505 

Maness H.L., Marcy G.W., Ford E.B., et al., 2007, PASP 1 19, 90 

Matsumura S., Takeda G, Rasio F.A., 2009, in "Transiting Planets", proc. IAU 

Symposium No 253, IAU, F. Pont, D. Sasselov, M. Holman (Eds.) 
Montagnier G, Bonfils X., Segransan D., et al., 2012, A&A, in preparation 
Pierens A., Nelson R.P, 2008, A&A 482, 333 

Pont E, Gilliland R.L., Knutson H.. Holman M.. Charbonneau D., 2009, 

MNRAS 393, L6 
Ribas I., Font-Ribera A., Beaulieu J.-R, 2008, ApJL 677, L59 
Saha P., Tremaine S., 1992, AJ 104, 1633 
Soderhjelm S., 1982, A&A 107, 54 
Tong Xiao, Zhou JiLin, 2009, Sci. Chin. Ser. G, 52, 4 
Torres G, 2007, ApJ 671, L65 
Trilling D.E., 2000, ApJ 537, L61 

Tsiganis K., Gomes R., Morbidelli A., Levison H E, 2005, Nature 435, 459 
Wisdom J., 2008, Icarus 193, 637 
Wu Y, Murray N., 2003, ApJ 589, 605 
Wu Y, Lithwick Y, 2011, ApJ 735, 109 
Zhang K, Hamilton D.P, 2008, Icarus 193, 267 



13 



